
%Edited by Katherine Hudson (November 21, 2014)
%Original Author: unknown

%To run this program, you must also have the following programs:
%exponentialfit.m, isodd.m, and SRfit.m

%Rename your data file data.txt. It's best to have a separate file for each
%data set. It's pretty awful. Blame the original author. I'm just polishing
%the turd here. 

%Program for Stress Relaxataion data, Poroelastic fits
clear all
batch = input('Do you want to batch process? Y/N [N]: ', 's');
if isempty(batch)
    batch = 'N';
end


folder_name = uigetdir;

cd(folder_name)




if strcmp(batch,'N')
    %select file to analyze
    [FileName,PathName,FilterIndex] = uigetfile([folder_name '/*.txt']);
    DATAfiles(1).name = FileName;
else
    DATAfiles = dir('*.txt');
    
end
nFiles = size(DATAfiles,1);


for kk = 1:nFiles
    FolderYN = exist(DATAfiles(kk).name,'dir');
    FolderName = DATAfiles(kk).name(1:end-4); %remove the extension
    if FolderYN == 0
        mkdir(FolderName); 
    end

sprintf('Processing file: %s',DATAfiles(kk).name)
%First step: User Define and Input Sample Info
diameter = input('Diameter of sample in mm:');
stepsize = input('Step displacement in mm(Enter 1):');

%% initialization of the code
%ini_disp = input('Initial Displacement in Wintest in m');
h0 = input('Height of sample in mm:');
nsteps = input('Number of steps:');



%Load data file
[h,data1] = hdrload(DATAfiles(kk).name); % data=[points, elapsedtime, scantime ,Displacement(mm), load(g)]
data = zeros(size(data1,1),5);
data(:,2) = data1(:,1);
data(:,4) = data1(:,2) - data1(1,2);
data(:,5) = data1(:,3);

ipt = findchangepts(data1(:,2),'MaxNumChanges',nsteps);

%save original data for Stressrelaxfit portion
dataStressfit = data;



%Parse data into time, stress, strain
time = data(:,2);
stress=data(:,5);%.*(.001*9.81*4)/(pi*(diameter/1000)^2); %(:,5)
step = [];%step(:,4)= -(data(:,4)./stepsize);
ini_disp = data(1,4);
% data(:,4)=data(:,4)-ini_disp;
strain=(data(:,4))./h0;     %(:,4)
ri = zeros(size(data1(:,2),1),1);%ri=abs(data(:,4)./stepsize);   %rough incrimenting (for chunking)  (:,4)

ri(1:ipt(1)-1) = 0;
for i=1:size(ipt,1)-1
ri(ipt(i):ipt(i+1)) = i ;
   
end
ri(ipt(end):end) = i+1;
data = [time, stress, strain, ri];

% (FB) this is already matlab default !(so why set it???)
format short


% repeat the last row of data for no reason and fix the index  (FB??)
data = [data; data(end,:)];
data(end,end) = nsteps + 1;


% round (:,4) to the nearest integer, throwing out rows not within .25
data = data(find(abs(round(data(:,4))-data(:,4)) < 0.25),:);
data(:,4) = round(data(:,4));


% find mean stress when ri=(:,4)=0
mean_stress = mean(data(find(data(:,4)==0),2));


% normalize stresses to mean stress
data(:,2) = data(:,2) - mean_stress;


% remove rows
data = data(find(data(1:end-1,2) < data(2:end,2)),:);

%% processing

%Change directory to Steps Folder....  that is the most inificient way to
%do it since it confuze matlab search path
PathSteps = [FolderName '/Steps'];
PathAnalysis = [FolderName '/Analysis'];
PathExist = exist(PathSteps,'dir');
if PathExist == 0
     mkdir(PathSteps);
     mkdir(PathAnalysis);
end


%Chunk data for each step

i=1;
index = zeros(nsteps,1);

for i = 1:nsteps

    I = find(data(:,4)==i);

    i;

    index(i) = I(1);
    index(i+1) = I(end)+1;

    chunkdata = data(I,1:3);

    filename = sprintf('%d.txt',i);
    fid = fopen([PathSteps '/' filename],'w');
    fprintf(fid,'%6.2f\t %12.8f\t %12.8f\n',chunkdata');
    status = fclose(fid);

    subplot(ceil(nsteps/2),2,i)
    plot(chunkdata(:,1),chunkdata(:,2),'.')
    xlabel('Time(s)')
    ylabel('Load (g)')
    title(sprintf('%d step',i));

end


%Ask user if the plots are ok

answer=input('Are you happy with the plots?  If so, enter 0.  Otherwise, enter the step number you would like to change\n');

while answer % if b is nonzero

    if answer > nsteps | answer < 1 ;

        answer = input('\n\nYou have entered a step plot that does not exist! \n\n Please try again.');

    else

        left = input ('Add how many points to the left?  (negative to subtract)\n');
        right = input ('Add how many points to the right?  (negative to subtract)\n');

        index(answer) = index(answer) - left;
        index(answer + 1) = index(answer + 1) + right;

        for i = answer - 1:answer + 1

            if i > 0 && i <= nsteps

                subplot(ceil(nsteps/2),2,i)
                I = index(i):index(i+1)-1;
                plot(data(I,1),data(I,2),'.')
                title(sprintf('%d step',i));

            end
        end
    end
    answer=input('Are you happy with the plots?  If so, enter 0.  Otherwise, enter the step number you would like to change\n');
end

disp('On to the Fit!');


saveas(gcf,[PathAnalysis '/graph of chunked data'],'fig')
saveas(gcf,[PathAnalysis '/graph of chunked data'],'pdf')




format compact

% mydir = pwd;
% addpath(pwd);
% cd(mydir);

dtop = dir;
sz = size(dtop,1);

%Fit each step

% cd steps;

for i = 1:nsteps

    name = sprintf([PathSteps '/%d.txt'],i);
%     dd=pwd;
    %disp(sprintf('%s:%s',dd,name));



    [model,real,values,fit]=SRfit(name);

    Stress(i,1:4)=model;         %[A,B,C,Strain]    A,B,C are Rawly's curve fit values-have to convert
    R2(i,1)=fit;
    numberplots=nsteps;

    if 1==isodd(nsteps)

        numberplots=nsteps+1;

    end

    subplot(numberplots/2,2,i)
    tit = sprintf('%d step',i);
    plot(real(:,1),real(:,2),'o',values(:,1),values(:,2),'rx-');
    legend('Data', 'Aexp(B*x)+C');





end



A=-1*Stress(:,1);

Tau=-(1./Stress(:,2));

B=Stress(:,3); %-Stress(:1)

stress2=[A,Tau,B,Stress(:,4),R2];



% cd ..

% cd Analysis



fid= fopen([PathAnalysis '/Parameters.txt'],'w');

fprintf(fid,'%6.2f\t  %12.8f\t %12.8f\t %12.8f\t %12.8f\n',stress2');

status=fclose(fid);

saveas(gcf,[PathAnalysis '/model fits graph'],'fig')
saveas(gcf,[PathAnalysis '/model fits graph'],'pdf')

% cd ..



%[h,info]=hdrload('info.txt');



% cd Analysis;

param=importdata([PathAnalysis '/Parameters.txt']);



tau=param(:,2);

eqstress=param(:,3);

inststress=(param(:,1)-param(:,3));

strain=param(:,4);



% graph the Eq Stress-Strain and Inst Stress-strain curve and ask user to select linear region

figure

subplot(2,1,1)

plotstrain = -1*strain;

ploteqstress = -1*eqstress;

plotinststress = -1*inststress;

plot(plotstrain,ploteqstress,'o');

subplot(2,1,2)

plot(plotstrain,inststress,'x');



ans=input('if all data looks good enter 0.  If not enter points you want to exclude (in row vector if more than one point \n[#,#,#]\n )');



while ans ~= 0

    tau(ans',:)=[];

    eqstress(ans',:)=[];

    inststress(ans',:)=[];

    strain(ans',:)=[];

    figure(2);

    subplot(2,1,1)

    plot(strain,eqstress,'o');

    subplot(2,1,2)

    plot(strain,inststress,'x');



    ans=input('If figure 2 is correct enter 0.\n if not re enter the points from figure 1 that you would like to remove.\n example[#,#,#]');

end



%plot(strain,eqstress,'o');

disp('Choose the linear region of the curve');

L = input('Lower bounds >>>  ');

U = input('Upper bounds >>>  ');

if L > U;

    disp('Error: Lower bound is not lower than Upper');

    disp('Exiting');

    pause(8);

    exit

end



figure

plot(dataStressfit(:,2),dataStressfit(:,5))

disp('Choose section for time constant analysis');

Lower = input('Lower time >>>  ');

Upper = input('Upper time >>>  ');

if Lower > Upper;

    disp('Error: Lower time is not lower than Upper time');

    disp('Exiting');

    pause(8);

    exit

end



% elimate unwanted rows

B = dataStressfit(dataStressfit(:,2)<=Upper,:);

B = B(B(:,2)>=Lower,:);



B(:,5) = -1*B(:,5);

B(:,2) = B(:,2)-Lower;



D = min(B);



B(:,5) = B(:,5)-(D(1,5)+.01);



%plot(B(:,2),B(:,5))



exponentialfit(B(:,2),B(:,5));



x = 0:(U-L);

y = ans(1).* exp(-ans(2) * x);

lambda = ans(2);



%

plot(B(:,2),B(:,5),'o',x,y,'-')





% Fit a linear regression to the range specified by user

Fstrain = strain(L:U,1);

Feqstress = eqstress(L:U,1);

p=polyfit(Fstrain,Feqstress,1);

f=polyval(p,Fstrain);

Finststress = inststress(L:U,1);

p2=polyfit(Fstrain,Finststress,1);

f2=polyval(p2,Fstrain);

clf;



% Get value for Ha, Inst and calculate permeability for each step



for i=1:length(tau)

    Ha(i)= (eqstress(i)- p(2))/(strain(i));

    Inst(i)=(inststress(i)-p2(2))/(strain(i));

    perm(i)=(((diameter/1000)/2)^2)./(pi^2*Ha(i)*tau(i));

end

perm=perm';



plotInst = -1*Inst;

plotFstrain = -1*Fstrain;

plotf = -1*f;



%Plot the curve with the fit and the permeability data

subplot(3,1,1);

plot(plotstrain,ploteqstress,'o',plotFstrain,plotf);

title('Eq Stress vs Strain');

%subplot(5,1,2);

%plot(plotstrain,perm,'kd',plotstrain,perm,'r-');

%title('Permeability vs Strain');

%subplot(5,1,3);

%plot(plotstrain,tau,'ks',plotstrain,tau,'b-');

%title('Tau vs Strain');

subplot(3,1,2);

plot(plotstrain,Ha, 'kp', plotstrain, Ha, 'g-');

title ('Equilibrium Modulus vs. strain');

subplot(3,1,3);

plot(plotstrain, plotInst, 'kp',plotstrain, plotInst);

title ('Instantaneous Modulus vs. strain');

saveas(gcf,[PathAnalysis '/HapermtauInst'],'fig');
saveas(gcf,[PathAnalysis '/HapermtauInst'],'pdf');


permfinal =(((diameter/1000)/2)^2)./(pi^2*p(1)*(1/lambda));



instant = -1*p2(1);

results=[param(:,1),param(:,2),param(:,3)];

fprintf('Equilibrium Modulus (Pa) %12.2f\n', p(1))

fprintf('Instanteous Modulus (Pa) %12.2f\n', instant)

fprintf('Permeability %17e\n', permfinal)



%constants=[strain,eqstress,perm];

fid= fopen([PathAnalysis '/Analysis.txt'],'w');

fprintf(fid,'%6.2f\t %12.8f\t %12.8f\t %12.8f\t  %12.8f\t %12.8f\n',results');

%fprintf(fid,'%6.2f\t %12.8f\t %12.8f\t %12.8f\t  %12.8f\t %12.8f\n',constants');%why does it only save 4 decimal places?

status=fclose(fid);

fid2=fopen([PathAnalysis '/Ha.txt'], 'w');

fprintf(fid, '%6.2f\t %12.8f\t %12.8f\t %12.8f\t  %12.8f\t %12.8f\n', p(1));

status=fclose(fid);

fid3=fopen([PathAnalysis '/Inst.txt'], 'w');

fprintf(fid, '%6.2f\t %12.8f\t %12.8f\t %12.8f\t  %12.8f\t %12.8f\n', instant);

status=fclose(fid);



% cd ..
close all
end